//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
gen c2 = cum_avg_daycare_36m_sum

summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

global inputs_std cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs c2

rename cum_avg_daycare_36m_sum c
rename parenting_ages13_std    p

// distribution of inputs and outputs
foreach var of varlist p c {
	foreach sib in 0 1 {
		bootstrap, strata(bwg site tg) reps($bootstraps): reg `var' tg if twin == `sib', robust 
		matrix b = e(b)
		matrix V = e(V)
		matrix  mc`var'_`sib' = b[1,2]
		matrix  md`var'_`sib' = b[1,1] 
		matrix se`var'_`sib' = sqrt(V[1,1])
	}
}

// labels
foreach var of varlist p c {
	foreach sib in 0 1 {
		foreach stat in mc md se {
			local `stat'`var'_`sib' = `stat'`var'_`sib'[1,1]
			local `stat'`var'_`sib' = string(``stat'`var'_`sib'',"%9.2f")
		}
	}
}

// only for childcare
foreach sib in 0 1 {
	summ   c2 if tg == 0 & twin == `sib'
	matrix sdcc_`sib'  = r(sd)
	local  sdcc_`sib'  = sdcc_`sib'[1,1]
	local  sdcc_`sib'  = string(`sdcc_`sib'',"%9.2f")
}

// parenting
cd $output
foreach sib in 0 1 {
	# delimit
	global boxp`sib'  text(.375 97.5
		 "{bf: Control Mean}"
		 "{bf: `mcp_`sib''}"
		 " "
		 "{bf: {&Delta}}" 
		 "{bf: `mdp_`sib'' (s.e. = `sep_`sib'')}"
	, place(c) box size(medsmall) just(c) margin(l+1 b+1 t+1 r+1) fcolor(none)); 
	# delimit cr

	#delimit
	twoway      (kdensity p if tg == 0 & twin == `sib', lwidth(dash) lcolor(gs0) lpattern(dash) lwidth(thick) yline(.45, lpattern(solid) lcolor(gs14)))
		    (kdensity p if tg == 1 & twin == `sib', lcolor(gs6)  lwidth(thick) yline(0, lpattern(solid) lcolor(gs14)) ),	  
			  

			legend(label(1 "Control") label(2 "Treatment") position(6) size(medsmall) order(1 2) rows(1) region(lcolor(gs0)))
			  xlabel(96[2]102, labsize(medium)        grid   glcolor(gs12) glpattern(solid)) 
			  ylabel(0[.15].45,   labsize(medium) angle(h) glcolor(gs12) glpattern(solid)) 
			  ytitle("Density", size(medium))
			  xtitle("Factor Score")
			  graphregion(color(white)) plotregion(color(white)) ${boxp`sib'};
	#delimit cr
	graph export parentingdist_s`sib'.eps, replace
}

// childcare
foreach sib in 0 1 {
	# delimit
	global boxc`sib'  text(.65 10.5
			 "{bf: Control Mean}"
			 "{bf: `mcc_`sib''}"
			 " "
			 "{bf: Control Standard Deviation}"
			 "{bf: `sdcc_`sib''}"
			 " "
			 "{bf: {&Delta}}" 
			 "{bf: `mdc_`sib'' (s.e. = `sec_`sib'')}"
	, place(c) box size(medsmall) just(c) margin(l+1 b+1 t+1 r+1) fcolor(none)); 
	# delimit cr

	#delimit
	twoway          (histogram c if tg == 1 & twin == `sib',  width(1)  fraction bcolor(gs6) gap(5) xline(12.3, lpattern(solid) lcolor(gs14)))
			(histogram c if tg == 0 & twin == `sib',  width(1)  fraction bfcolor(none) blcolor(gs0)  blwidth(medthick)  gap(5)),	  
			  

			legend(label(2 "Control") label(1 "Treatment") position(6) size(medsmall) order(2 1) rows(1) region(lcolor(gs0)))
			  xlabel(5.3[1]12.3, labsize(medium)        grid   glcolor(gs12) glpattern(solid)) 
			  ylabel(0[.30].90,   labsize(medium) angle(h) glcolor(gs12) glpattern(solid)) 
			  ytitle("Fraction", size(medium))
			  xtitle("Standardized Average Hours of Childcare per Week")
			  graphregion(color(white)) plotregion(color(white)) ${boxc`sib'};
	#delimit cr
	graph export childcaredist_s`sib'.eps, replace
}
